function [beta,R2] = reg_Y_NxYzXe_on_X_Nxz(Y_Nxz,coeffY_Y,coeffY_X,coeffY_e,grid_Y,...
    inflation_process,X_Nxz,prob_NxYz)
% note: Y is specialized to inputs related to nominal bonds
% Y = Y_Nxz + coeffY_Y*Y + coeffY_X*X + coeffY_e*e
% (X,e) is independent from (N,x,Y,z)

    mean_X = inflation_process.mu_pi/(1 - inflation_process.rho_pi);
    var_X = inflation_process.sigma_pi^2*(1 ...
        + 2*inflation_process.rho_pi*inflation_process.theta_pi ...
        + inflation_process.theta_pi^2)/(1 - inflation_process.rho_pi^2);
    var_e = inflation_process.sigma_pi^2;
    cov_Xe = inflation_process.sigma_pi^2;

    [NN, Nx, NY, Nz] = size(prob_NxYz);

    Y_NxYz = permute(repmat(Y_Nxz,[1 1 1 NY]),[1 2 4 3]) ...
        + coeffY_Y*permute(repmat(grid_Y(:),[1 NN Nx Nz]),[2 3 1 4]);
    mean_Y_NxYz = sum(prob_NxYz.*Y_NxYz,'all');
    var_Y_NxYz = sum(prob_NxYz.*((Y_NxYz - mean_Y_NxYz).^2),'all');
    mean_lhs = mean_Y_NxYz + coeffY_X*mean_X;
    var_lhs = var_Y_NxYz + coeffY_X^2*var_X + coeffY_e^2*var_e ...
        + 2*coeffY_X*coeffY_e*cov_Xe;

    prob_Nxz = squeeze(sum(prob_NxYz,3));
    mean_rhs = sum(X_Nxz.*prob_Nxz,'all');
    var_rhs = sum(((X_Nxz - mean_rhs).^2).*prob_Nxz,'all');

    E_lhs_rhs = sum(prob_NxYz.*Y_NxYz.*permute(repmat(X_Nxz,[1 1 1 NY]),[1 2 4 3]),'all') ...
        + coeffY_X*mean_X*mean_rhs;
    cov_lhs_rhs = E_lhs_rhs - mean_lhs*mean_rhs;

    beta = cov_lhs_rhs/var_rhs;
    R2 = beta^2*var_rhs/var_lhs;

end